Language/OS - Multiplatform Resource Library

home *** CD-ROM | disk | FTP | other *** search

/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / cpp_libs / newmat03.lha / newmat03 / newmat.txt < prev next >

Wrap

Text File | 1993-08-08 | 53KB | 1,500 lines

//$$ newmat.txt Documentation file Documentation for newmat03, an experimental matrix package in C++. ================================================================== MATRIX PACKAGE 25 November, 1991 Copyright (C) 1991: R B Davies and DSIR Permission is granted to use but not to sell. Contents ======== General description Is this the package you need? Changes Where you can get a copy of this package Compiler performance Example Detailed documentation Customising Constructors Elements of matrices Matrix copy Unary operators Binary operators Combination of a matrix and scalar Scalar functions of matrices Submatrix operations Change dimensions Change type Multiple matrix solve Memory management Output Accessing matrices of unspecified type Cholesky decomposition Householder triangularisation Singular Value Decomposition Eigenvalues Sorting Fast Fourier Transform Interface to Numerical Recipes in C List of files Notes on the design of the package What this is package for What size of matrices? Allow matrix expressions? Which matrix types? What element types? Naming convention Row and Column index ranges Structure of matrix objects Data storage - one block or several Data storage - by row or by column or other Storage of symmetric matrices Element access - method and checking Use iterators? Memory management - reference counting or status variable? Evaluation of expressions - use two stage method? How to overcome an explosion in number of operations Using const A calculus of matrix types Error handling Band and sparse matrices Problem report form --------------------------------------------------------------------------- General description =================== The package is intented for scientists and engineers who need to manipulate a variety of types of matrices using standard matrix operations. Emphasis is on the kind of operations needed in statistical calculations such as least squares, linear equation solve and eigenvalues. It supports matrix types Matrix (rectangular matrix) nricMatrix (variant of rectangular matrix) UpperTriangularMatrix LowerTriangularMatrix DiagonalMatrix SymmetricMatrix RowVector (derived from Matrix) ColumnVector (derived from Matrix). Only one element type (float or double) is supported. The package includes the operations *, +, -, inverse, transpose, conversion between types, submatrix, determinant, Cholesky decomposition, Householder triangularisation, singular value decomposition, eigenvalues of a symmetric matrix, sorting, fast fourier transform, printing and an interface with "Numerical Recipes in C". It is intended for matrices in the range 4 x 4 to about 90 x 90 (125 x 125 for triangular matrices). The upper limit is imposed by the maximum number of elements that can be contained in a single array (8192 doubles in some machines). A two-stage approach to evaluating matrix expressions is used to improve efficiency and reduce use of temporary storage. The package is designed for version 2 of C++. It works with Turbo C++, Borland C++, Glockenspiel C++ (2.00a) on a PC and AT&T C++ (2.0) and Gnu C++ on a Sun. It works with some problems with Zortech C++ (version 2). --------------------------------------------------------------------------- Is this the package you need? ============================= Do you 1. need matrix operators such as * and + defined as operators so you can write things like X = A * (B + C); 2. need a variety of types of matrices 3. need only one element type (float or double) 4. work with matrices in the range 4x4 to 90x90 5. tolerate a large and complex package Then maybe this is the right package for you. If you don't need (1) then there may be better options. Likewise if you don't need (2) there may be better options. If you require "not (5)" then this is not the package for you. If you need (2) and "not (3)" and have some spare money, then maybe you should look at M++ from Dyad or the Rogue Wave matrix package. --------------------------------------------------------------------------- Changes ======= Newmat03 - November 1991: Col and Cols become Column and Columns. Added Sort, SVD, Jacobi, Eigenvalues, FFT, real conversion of 1x1 matrix, "Numerical Recipes in C" interface, output operations, various scalar functions. Improved return from functions. Reorganised setting options in "include.hxx". Newmat02 - July 1991: Version with matrix row/column operations and numerous additional functions. Matrix - October 1990: Early version of package. --------------------------------------------------------------------------- How to get a copy of this package ================================= I am putting copies on Compuserve (Borland library, zip format), SIMTEL20 (MsDos library, zip format), comp.sources.misc on Internet (shar format), and on the MsDos program library at Victoria University, Wellington. --------------------------------------------------------------------------- Compiler performance ==================== I have tested this package on a number of compilers. Here are the levels of success with this package. In most cases I have chosen code that works under all the compilers I have access to, but I have had to include some specific work-arounds for some compilers. For the MsDos versions, I am using a 386/387sx computer running MsDos 5, except that Turbo is on an old XT. The unix versions are on a Sun Sparc station. A series of #defines at the beginning of "include.hxx" customises the package for the compiler you are using. Turbo, Borland and Zortech are recognised automatically, otherwise you have to set the appropriate #define statement. The compilers are looking a bit old now. I do intend to test the package against newer versions as they become available. Borland C++ 2.0: Recently, this has been my main development platform, so naturally almost everything works with this compiler. The library manager "tlib" fails but you can use "zorlib" from Zortech instead. Sometimes Borland crashes during a compiler or mis-compiles. You just have to reboot and continue the compile. Turbo C++ (? version): Almost works OK. My rather elderly version does show a problem. Probably not worth tracking down - buy a newer version. Haven't tried the linker. Zortech C++ 2.12: "const" doesn't work correctly with this compiler, so the package skips all of the statements Zortech can't handle. If you are using a later version of Zortech you could probably re-activate these statements. Zortech leaves rubbish on the heap. I don't know whether this is my programming error or a Zortech error. It works better when one doesn't optimise but there still are problems. The nric routines don't work. Zortech does not support IO manipulators. Glockenspiel C++ (2.00a for MsDos loading into Microsoft C 5.1): I haven't tested the latest version of my package with Glockenspiel. I had to #define the matrix names to shorter names to avoid ambiguities and had quite a bit of difficulty stopping the compiles from running out of space and not exceeding Microsoft's block nesting limit. A couple of my test statements produced statements too complex for Microsoft, but basically the package worked. This was my original development platform and I still use .cxx and .hxx as my file name extensions. However, Glockenspiel is no longer competitive for MsDos and I am not updating my copy of the compiler. Sun AT&T C++ 2.00: This works fine. Except aggregates are not supported. Gnu G++ 1.37.1: This mostly works. You don't seem to be able to use expressions like Matrix(X*Y) in the middle of an expression and (Matrix)(X*Y) is unreliable. Gnu does not support IO manipulators. Gnu leaves rubbish on the heap. This is from output statements and not my package and may not be an error. The previous version of the package did not work under Gnu 1.37 or 1.39. --------------------------------------------------------------------------- Example ======= An example is given in example.cxx . This gives a simple linear regression example using four different algorithms. The correct output is given in example.txt. The program carries out a check that no memory is left allocated on the heap when it terminates. The file example.dep contains a dependency list for compiling example.cxx . You will need to add the compile and link commands. example.dep lists all the files in the package so you can adapt for other projects. Don't forget to remove references to newmat9.cxx if you are using a compiler that does not support the standard io manipulators. --------------------------------------------------------------------------- Detailed Documentation ====================== Copyright (C) 1989,1990,1991: R B Davies and DSIR Permission is granted to use but not to sell. -------------------------------------------------------------- | Please understand that this is a test version; there may | | still be bugs and errors. Use at your own risk. Neither I | | nor DSIR take any responsibility for any errors or omissions | | in this package or for any misfortune that may befall you or | | others as a result of its use. | -------------------------------------------------------------- Please report bugs to me at robert@am.dsir.govt.nz or Compuserve 72777,656 When reporting a bug please tell me which C++ compiler you are using (if known), and what version. Also give me details of your computer (if known). Tell me where you downloaded your version of my package from and its version number (eg newmat03 or newmat04). (There may be very minor differences between versions at different sites). Note any changes you have made to my code. If at all possible give me a piece of code illustrating the bug. Please do report bugs to me. The matrix inverse routine and the sort routines are adapted from "Numerical Recipes in C" by Press, Flannery, Teukolsky, Vetterling, published by the Cambridge University Press. Other code is adapted from routines in "Handbook for Automatic Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published by Springer Verlag. Customising ----------- I use .hxx as the suffix of definition files and .cxx as the suffix of C++ source files. This does not cause any problems with the compilers I use except that Borland and Turbo need to be told to accept any suffix as meaning a C++ file rather than a C file. Use the large model when you are using a PC. Do not "outline" inline functions. Each file accessing the matrix package needs to have file newmat.hxx #included at the beginning. Files using matrix applications (Cholesky decomposition, Householder triangularisation) need newmatap.hxx instead (or as well). If you need the output functions you will also need newmatio.hxx. Glockenspiel also needs to have include.hxx #included before newmat.hxx. The file include.hxx sets the options for the compiler. If you are using a compiler different from one I have worked with you may have to set up a new section in include.hxx appropriate for your compiler. Borland, Turbo and Zortech are recognised automatically. If you using Glockenspiel on a PC, AT&T, or Gnu C++ activate the appropriate statement at the beginning of include.hxx. Activate the appropriate statement to make the element type float or double. The file (newmat9.cxx) containing the output routines can be used only with libraries that support the AT&T input/output routines including manipulators. It cannot be used with Zortech or Gnu. Constructors ------------ To construct an m x n matrix, A, (m and n are integers) use Matrix A(m,n); The UpperTriangularMatrix, LowerTriangularMatrix, SymmetricMatrix and DiagonalMatrix types are symmetric. To construct an n x n matrix use, for example UpperTriangularMatrix U(n); Likewise the RowVector and ColumnVector types take just one argument in their constructors: RowVector RV(n); You can also construct vectors and matrices without specifying the dimension. For example Matrix A; In this case the dimension must be set by an assignment statement or a re-dimension statement. You can also use a constructor to set a matrix equal to another matrix or matrix expression. Matrix A = U; Matrix A = U * L; Only conversions that don't lose information are supported - eg you cannot convert an upper triangular matrix into a diagonal matrix using =. Elements of matrices -------------------- Elements are accessed by expressions of the form A(i,j) where i and j run from 1 to the appropriate dimension. Access elements of vectors with just one argument. Diagonal matrices can accept one or two subscripts. This is different from the earlier version of the package in which the subscripts ran from 0 to one less than the appropriate dimension. Use A.element(i,j) if you want this earlier convention. Matrix copy ----------- The operator = is used for copying matrices, converting matrices, or evaluating expressions. For example A = B; A = L; A = L * U; Only conversions that don't lose information are supported. The dimensions of the matrix on the left hand side are adjusted to those of the matrix or expression on the right hand side. Elements on the right hand side which are not present on the left hand side are set to zero. The operator << can be used in place of = where it is permissible for information to be lost. For example SymmetricMatrix S; Matrix A; ...... S << A.t() * A; is acceptable whereas S = A.t() * A; // error will cause a runtime error since the package does not (yet) recognise A.t()*A as symmetric. Note that you can not use << with constructors. For example SymmetricMatrix S << A.t() * A; // error does not work. A third copy routine is used in a similar role to =. Use A.Inject(D); to copy the elements of D to the corresponding elements of A but leave the elements of A unchanged if there is no corresponding element of D (the = operator would set them to 0). This is useful, for example, for setting the diagonal elements of a matrix without disturbing the rest of the matrix. Unlike = and <<, Inject does not reset the dimensions of A, which must match those of D. Both << and Inject can be used with submatrix expressions on the left hand side. See the section on submatrices. To set the elements of a matrix to a scalar use operator = real r; Matrix A(m,n); ...... Matrix A(m,n); A = r; You can load the elements of a matrix from an array: Matrix A(3,2); real a[] = { 11,12,21,22,31,33 }; A << a; This construction cannot check that the numbers of elements match correctly. This version of << can be used with submatrices on the left hand side. Unary operators --------------- The package supports unary operations change sign of elements -A transpose A.t() inverse (of square matrix A) A.i() Binary operations ----------------- The package supports binary operations matrix addition A+B matrix subtraction A-B matrix multiplication A*B equation solve (square matrix A) A.i()*B In the last case the inverse is not calculated. Notes: If you are doing repeated multiplication. For example A*B*C, use brackets to force the order to minimize the number of operations. If C is a column vector and A is not a vector, then it will usually reduce the number of operations to use A*(B*C) . The package does not recognise B*A.i() as an equation solve. It is probably better to use (A.t().i()*B.t()).t() . Combination of a matrix and scalar ---------------------------------- The following expression multiplies the elements of a matrix A by a scalar f: A * f; Likewise one can divide the elements of a matrix A by a scalar f: A / f; The expressions A + f and A - f add or subtract a rectangular matrix of the same dimension as A with elements equal to f to or from the matrix A. In each case the matrix must be the first term in the expression. Expressions such f + A or f * A are not recognised. Scalar functions of matrices ---------------------------- int m = A.Nrows(); // number of rows int n = A.Ncols(); // number of columns real ssq = A.SumSquare(); // sum of squares of elements real sav = A.SumAbsoluteValue(); // sum of absolute values real mav = A.MaximumAbsoluteValue(); // maximum of absolute values real norm = A.Norm1(); // maximum of sum of absolute values of elements of a column real norm = A.NormInfinity(); // maximum of sum of absolute values of elements of a row real t = A.Trace(); // trace LogandSign ld = A.LogDeterminant(); // log of determinant BOOL z = A.IsZero(); // test all elements zero MatrixType mt = A.Type(); // type of matrix real* s = Store(); // pointer to array of elements int l = Storage(); // length of array of elements A.LogDeterminant() returns a value of type LogandSign. If ld is of type LogAndSign use ld.Value() to get the value of the determinant ld.Sign() to get the sign of the determinant (values 1, 0, -1) ld.LogValue() to get the log of the absolute value. A.IsZero() returns BOOL value TRUE if the matrix A has all elements equal to 0.0. MatrixType mt = A.Type() returns the type of a matrix. Use (char*)mt to get a string (UT, LT, Rect, Sym, Diag, RowV, ColV, Crout) showing the type. SumSquare(A), SumAbsoluteValue(A), MaximumAbsoluteValue(A), Trace(A), LogDeterminant(A), Norm1(A), NormInfinity(A) can be used in place of A.SumSquare(), A.SumAbsoluteValue(), A.MaximumAbsoluteValue(), A.Trace(), A.LogDeterminant(), A.Norm1(), A,NormInfinity(). Submatrix operations -------------------- A.SubMatrix(fr,lr,fc,lc) This selects a submatrix from A. the arguments fr,lr,fc,lc are the first row, last row, first column, last column of the submatrix with the numbering beginning at 1. This may be used in any matrix expression or on the left hand side of << or Inject. Inject does not check no information loss in this case. You can also use the construction real c; .... A.SubMatrix(fr,lr,fc,lc) << c; to set a submatrix equal to a constant. The follwing are variants of SubMatrix: A.SymSubMatrix(f,l) // This assumes fr=fc and lr=lc. A.Rows(f,l) // select rows A.Row(f) // select single row A.Columns(f,l) // select columns A.Column(f) // select single column In each case f and l mean the first and last row or column to be selected (starting at 1). If SubMatrix or its variant occurs on the right hand side of an = or << or within an expression its type is as follows A.Submatrix(fr,lr,fc,lc): If A is RowVector or ColumnVector then same type otherwise type Matrix A.SymSubMatrix(f,l): Same type as A A.Rows(f,l): Type Matrix A.Row(f): Type RowVector A.Columns(f,l): Type Matrix A.Column(f): Type ColumnVector If SubMatrix or its variant appears on the left hand side of << , think of its type being Matrix. Thus L.Row(1) where L is LowerTriangularMatrix expects L.Ncols() elements even though it will use only one of them. Change dimensions ----------------- The following operations change the dimensions of a matrix. The values of the elements are lost. A.ReDimension(nrows,ncols); // for type Matrix or nricMatrix A.ReDimension(n); // for all other types Change type ----------- The following functions interpret the elements of a matrix (stored row by row) to be a vector or matrix of a different type. Actual copying is usually avoided where these occur as part of a more complicated expression. A.CopyToRow() A.CopyToColumn() A.CopyToDiagonal() A.CopyToMatrix(nrows,ncols) A.c() real(A) The form .c() is used in matrix expressions when A is of a const type. The expression real(A) is used to convert a 1 x 1 matrix to a scalar. Multiple matrix solve --------------------- If A is a square or symmetric matrix use CroutMatrix X = A; // carries out LU decomposition Matrix AP = X.i()*P; Matrix AQ = X.i()*Q; LogAndSign ld = X.LogDeterminant(); rather than Matrix AP = A.i()*P; Matrix AQ = A.i()*Q; LogAndSign ld = A.LogDeterminant(); since each operation will repeat the LU decompostion. Memory management ----------------- The package does not support delayed copy. Several strategies are required to prevent unnecessary matrix copies. Where a matrix is called as a function argument use a constant reference. For example YourFunction(const Matrix& A) rather than YourFunction(Matrix A) Constant matrices cannot be used in matrix expressions so if you wish to use A in an expression within this function use A.c() rather than A. Skip the rest of this section on your first reading. A second place where it is desirable to avoid unnecessary copies is when a function is returning a matrix. Matrices can be returned from a function with the return command as you would expect. However these may incur one and possibly two copyings of the matrix. To avoid this use the following instructions. Make your function of type ReturnMatrix . Then precede the return statement with a Release statement (or a ReleaseAndDelete statement if the matrix was created with new). For example ReturnMatrix MakeAMatrix() { Matrix A; ...... A.Release(); return A; } or ReturnMatrix MakeAMatrix() { Matrix* m = new Matrix; ...... m->ReleaseAndDelete(); return *m; } Note that .c() cannot be applied to a matrix following application of .Release() or ->ReleaseAndDelete() . --------------------------------------------------------------------- | Do not forget to make the function of type ReturnMatrix; otherwise | | incomprehensible run-time errors will occur with some compilers. | --------------------------------------------------------------------- You can also use .Release() or ->ReleaseAndDelete() to allow a matrix expression to recycle space. Suppose you call A.Release(); just before A is used just once in an expression. Then the memory used by A is either returned to the system or reused in the expression. In either case, A's memory is destroyed. This procedure can be used to imporve efficiency and reduce the use of memory. Use ->ReleaseAndDelete for matrices created by new if you want to completely delete the matrix after it is accessed. Output ------ To print a matrix use an expression like Matrix A; ...... cout << setw(10) << setprecision(5) << A; This will work only with systems that support the AT&T input/output routines including manipulators. Accessing matrices of unspecified type -------------------------------------- Skip this section on your first reading. Suppose you wish to write a function which accesses a matrix of unknown type including expressions (eg A*B). Then use a layout similar to the following: void YourFunction(BaseMatrix& X) { GeneralMatrix* gm = X.Evaluate(); // evaluate an expression // if necessary ........ // operations on *gm gm->tDelete(); // delete *gm if a temporary } See, as an example, the definitions of operator<< in newmat9.cxx. Cholesky decomposition ---------------------- Suppose S is symmetric and positive definite. Then there exists a unique lower triangular matrix L such that L * L.t() = S. To calculate this use SymmetricMatrix S; ...... LowerTriangularMatrix L = Cholesky(S); Householder triangularisation ----------------------------- Start with matrix / X 0 \ s \ Y 0 / t n s The Householder triangularisation post multiplies by an orthogonal matrix Q such that the matrix becomes / 0 L \ s \ Z M / t n s where L is lower triangular. Note that X is the transpose of the matrix sometimes considered in this context. This is good for solving least squares problems: choose b (matrix or row vector) to minimize the sum of the squares of the elements of Y - b*X Then choose b = M * L.i(); Two routines are provided: HHDecompose(X, L); replaces X by orthogonal columns and forms L. HHDecompose(X, Y, M); uses X from the first routine, replaces Y by Z and forms M. Singular Value Decomposition ---------------------------- The singular value decomposition of an m x n matrix A ( where m >= n) is a decomposition A = U * D * V.t() where U is m x n with U.t() * U equalling the identity, D is an n x n diagonal matrix and V is an n x n orthogonal matrix. Singular value decompositions are useful for understanding the structure of ill-conditioned matrices, solving least squares problems, and for finding the eigenvalues of A.t() * A. To calculate the singular value decomposition of A (with m >= n) use one of SVD(A, D, U, V); // U (= A is OK) SVD(A, D); SVD(A, D, U); // U (= A is OK) SVD(A, D, U, FALSE); // U (can = A) for workspace only SVD(A, D, U, V, FALSE); // U (can = A) for workspace only The values of A are not changed unless A is also inserted as the third argument. Eigenvalues ----------- An eigenvalue decomposition of a symmetric matrix A is a decomposition A = V * D * V.t() where V is an orthogonal matrix and D is a diagonal matrix. Eigenvalue analyses are used in a wide variety of engineering, statistical and other mathematical analyses. The package includes two algorithms: Jacobi and Householder. The first is extremely reliable but much slower than the second. The code is adapted from routines in "Handbook for Automatic Computation, Vol II, Linear Algebra" by Wilkinson and Reinsch, published by Springer Verlag. Jacobi(A,D,S,V); // A, S symmetric is for workspace, // S = A is OK Jacobi(A,D); // A symmetric Jacobi(A,D,S); // A, S symmetric is for workspace, // S = A is OK Jacobi(A,D,V); // A symmetric EigenValues(A,D); // A symmetric EigenValues(A,D,S); // A, S symmetric is for back // transforming, S = A is OK EigenValues(A,D,V); // A symmetric Sorting ------- To sort the values in a matrix or vector, A, (in general this operation makes sense only for vectors and diagonal matrices) use SortAscending(A); or SortDescending(A); Fast Fourier Transform ---------------------- FFT(CV1, CV2, CV3, CV4); // CV3=CV1 and CV4=CV2 is OK where CV1, CV2, CV3, CV4 are column vectors. CV1 and CV2 are the real and imaginary input vectors; CV3 and CV4 are the real and imaginary output vectors. The lengths of CV1 and CV2 must be equal and should be the product of numbers less than about 10 for fast execution. Interface to Numerical Recipes in C ----------------------------------- This package can be used with the vectors and matrices defined in "Numerical Recipes in C". You need to edit the routines in Numerical Recipes so that the elements are of the same type as used in this package. Eg replace float by double, vector by dvector and matrix by dmatrix, etc. You will also need to edit the function definitions to use the version acceptable to your compiler. Then enclose the code from Numerical Recipes in extern "C" { ... }. You will also need to include the matrix and vector utility routines. Then any vector in Numerical Recipes with subscripts starting from 1 in a function call can be replaced by a RowVector, ColumnVector or DiagonalMatrix in the present package. Similarly any matrix with subscripts starting from 1 can be replaced by an nricMatrix in the present package. The class nricMatrix is derived from Matrix and can be used in place of Matrix. Numerical Recipes cannot change the dimensions of a matrix or vector. So matrices or vectors must be correctly dimensioned before a Numerical Recipes routine is called. For example SymmetricMatrix B(44); ..... // load values into B nricMatrix BX = B; // copy values to an nricMatrix DiagonalMatrix D(44); // Matrices for output nricMatrix V(44,44); // correctly dimensioned int nrot; jacobi(BX,44,D,V,&nrot); // jacobi from NRIC cout << D; // print eigenvalues --------------------------------------------------------------------------- List of files ============= README readme file NEWMAT TXT documentation file NEWMAT LIS list of files CLASS TXT list of classes BOOLEAN HXX boolean class definition CONTROLW HXX control word definition file INCLUDE HXX details of include files and options NEWMAT HXX main matrix clss definition file NEWMATAP HXX applications definition file NEWMATIO HXX input/output definition file NEWMATRC HXX row/column functions definition files NEWMATRM HXX rectangular matrix access definition files PRECISIO HXX numerical precision constants CHOLESKY CXX Cholesky decomposition EVALUE CXX eigenvalues and eigenvector calculation FFT CXX fast Fourier transform HHOLDER CXX Householder triangularisation JACOBI CXX eigenvalues by the Jacobi method NEWMAT1 CXX type manipulation routines NEWMAT2 CXX row and column manipulation functions NEWMAT3 CXX row and column access functions NEWMAT4 CXX constructors, redimension, utilities NEWMAT5 CXX transpose, evaluate, matrix functions NEWMAT6 CXX operators, element access NEWMAT7 CXX invert, solve, binary operations NEWMAT8 CXX LU decomposition, scalar functions NEWMAT9 CXX output routines NEWMATRM CXX rectangular matrix access functions SORT CXX sorting functions SUBMAT CXX submatrix functions SVD CXX singular value decomposition EXAMPLE CXX example of use of package EXAMPLE TXT output from example EXAMPLE DEP dependency file for example --------------------------------------------------------------------------- Notes on the design of the package ================================== Copyright (C) 1991: R B Davies and DSIR Please treat this as an academic publication. You can use the ideas but please acknowledge in any publication. In this section, I describe some of the ideas behind this package, some of the decisions that I needed to make and give some details about the way it works. You don't need to read this section in order to use the package. I don't think it is obvious what is the best way of going about structuring a matrix package. And I don't think you can figure this out with "thought" experiments. Different people have to try out different approaches. And someone else may have to figure out which is best. Or, more likely, the ultimate packages will lift some ideas from each of a variety of trial packages. So, I don't claim my package is an "ultimate" package, but simply a trial of a number of ideas. But I do hope it will meet the immediate requirements of some people who need to carry out matrix calculations. What this is package for ------------------------ The package is used for the manipulation of matrices, including the standard operations such as multiplication as understood by numerical analysts, engineers and mathematicians. A matrix is a two dimensional array of numbers. However, very special operations such as matrix multiplication are defined specifically for matrices. This means that a "matrix" package tends to be different from a general "array" package. I see a matrix package as providing the following 1. Evaluation of matrix expressions in a form familiar to scientists and engineers. For example A = B * (C + D.t()); 2. Access to the elements of a matrix; 3. Access to submatrices; 4. Common elementary matrix functions such as determinant and trace; 5. A platform for developing advanced matrix functions such as eigen- value analysis; 6. Good efficiency and storage management; 7. Graceful exit from errors (I don't provide this yet). It may also provide 8. A variety of types of elements (eg real and complex); 9. A variety of types of matrices (eg rectangular and symmetric). In contrast an array package should provide 1'. Arrays can be copied with a single instruction; may have some arithmetic operations, say + - and scalar + - * /, it may provide matrix multiplication as a function; 2'. High speed access to elements directly and with iterators; 3'. Access to subarrays and rows (and columns?); 6'. Good efficiency and storage management; 7'. Graceful exit from errors; 8'. A variety of types of elements (eg real and complex); 9'. One, two and three dimensional arrays, at least, with starting points of the indices set by user; may have arrays with ragged rows. It may be possible to amalgamate these two sets of requirements to some extent. However my package is definitely oriented towards the first set. Even within the bounds set by the first set of requirements there is a substantial opportunity for variation between what different matrix packages might provide. It is not possible to build a matrix package that will meet everyone's requirements. In many cases if you put in one facility, you impose overheads on everyone using the package. This both is storage required for the program and in efficiency. Likewise a package that is optimised towards handling large matrices is likely to become less efficient for very small matrices where the administration time for the matrix may become significant compared with the time to carry out the operations. It is better to provide a variety of packages (hopefully compatible) so that most users can find one that meets their requirements. This package is intended to be one of these packages; but not all of them. Since my background is in statistical methods, this package is oriented towards the kinds things you need for statistical analyses. What size of matrices? ---------------------- A matrix package may target small matrices (say 3 x 3), or medium sized matrices, or very large matrices. A package targeting very small matrices will seek to minimise administration. A package for medium sized or very large matrices can spend more time on administration in order to conserve space or optimise the evaluation of expressions. A package for very large matrices will need to pay special attention to storage and numerical properties. This package is designed for medium sized matrices. This means it is worth introducing some optimisations, but I don't have to worry about the 64K limit that some compilers impose. Allow matrix expressions? ------------------------- I want to be able to write matrix expressions the way I would on paper. So if I want to multiply two matrices and then add the transpose of a third one I can write something like X = A * B + C.t(); I want this expression to be evaluated with close to the same efficiency as a hand-coded version. This is not so much of a problem with expressions including a multiply since the multiply will dominate the time. However, it is not so easy to achieve with expressions with just + and - . A second requirement is that temporary matrices generated during the evaluation of an expression are destroyed as quickly as possible. A desirable feature is that a certain amount of "intelligence" be displayed in the evaluation of an expression. For example, in the expression X = A.i() * B; where i() denotes inverse, it would be desirable if the inverse wasn't explicitly calculated. Which matrix types? ------------------- As well as the usual rectangular matrices, matrices occuring repeatedly in numerical calculations are upper and lower triangular matrices, symmetric matrices and diagonal matrices. This is particularly the case in calculations involving least squares and eigenvalue calculations. So as a first stage these were the types I decided to include. It is also necessary to have types row vector and column vector. In a "matrix" package, in contrast to an "array" package, it is necessary to have both these types since they behave differently in matrix expressions. The vector types can be derived for the rectangular matrix type, so having them does not greatly increase the complexity of the package. The problem with having several matrix types is the number of versions of the binary operators one needs. If one has 5 distinct matrix types then a simple package will need 25 versions of each of the binary operators. In fact, we can evade this problem, but at the cost of some complexity. What element types? ------------------- Ideally we would allow element types double, float, complex and int, at least. It would be reasonably easy, using templates or equivalent, to provide a package which could handle a variety of element types. However, as soon as one starts implementing the binary operators between matrices with different element types, again one gets an explosion in the number of operations one needs to consider. Hence I decided to implement only one element type. But the user can decide whether this is float or double. The package assumes elements are of type real. The user typedefs to float or double. In retrospect, it would not be too hard to include matrices with complex matrix type. It might also be worth including symmetric and triangular matrices with extra precision elements (double or long double) to be used for storage only and with a minimum of operations defined. These would be used for accumulating the results of sums of squares and product matrices or multistage Householder triangularisations. Naming convention ----------------- How are classes and public member functions to be named? As a general rule I have spelt identifiers out in full with individual words being capitalised. For example "UpperTriangularMatrix". If you don't like this you can #define or typedef shorter names. This convention means you can select an abbreviation scheme that makes sense to you. The convention causes problems for Glockenspiel C++ on a PC feeding into Microsoft C. The names Glockenspiel generates exceed the the 32 characters recognised by Microsoft C and ambiguities result. So it is necessary to #define shorter names. Exceptions to the general rule are the functions for transpose and inverse. To make matrix expressions more like the corresponding mathematical formulae, I have used the single letter abbreviations, t() and i() . Row and Column index ranges --------------------------- In mathematical work matrix subscripts usually start at one. In C, array subscripts start at zero. In Fortran, they start at one. Possibilities for this package were to make them start at 0 or 1 or be arbitrary. Alternatively one could specify an "index set" for indexing the rows and columns of a matrix. One would be able to add or multiply matrices only if the appropriate row and column index sets were identical. In fact, I adopted the simpler convention of making the rows and columns of a matrix be indexed by an integer starting at one, following the traditional convention. In an earlier version of the package I had them starting at zero, but even I was getting mixed up when trying to use this earlier package. So I reverted to the more usual notation. Structure of matrix objects --------------------------- Each matrix object contains the basic information such as the number of rows and columns and a status variable plus a pointer to the data array which is on the heap. Data storage - one block or several ----------------------------------- In this package the elements of the matrix are stored as a single array. Alternatives would have been to store each row as a separate array or a set of adjacent rows as a separate array. The present solution simplifies the program but limits the size of matrices in systems that have a 64k byte (or other) limit on the size of arrays. The large arrays may also cause problems for memory management in smaller machines. Data storage - by row or by column or other ------------------------------------------- In Fortran two dimensional arrays are stored by column. In most other systems they are stored by row. I have followed this later convention. This makes it easier to interface with other packages written in C but harder to interface with those written in Fortran. An alternative would be to store the elements by mid-sized rectangular blocks. This might impose less strain on memory management when one needs to access both rows and columns. Storage of symmetric matrices ----------------------------- Symmetric matrices are stored as lower triangular matrices. The decision was pretty arbitrary, but it does slightly simplify the Cholesky decomposition program. Element access - method and checking ------------------------------------ We want to be able to use the notation A(i,j) to specify the (i,j)-th element of a matrix. This is the way mathematicians expect to address the elements of matrices. I didn't even consider using the totally alien notation A[i][j]. There are two ways of working out the address of A(i,j). One is using a "dope" vector which contains the first address of each row. This is how C works when you use A[i][j]. Alternatively you can calculate the address using the formula appropriate for the structure of A. I use this second approach. It is probably slower, but saves worrying about an extra bit of storage. The other question is whether to check for i and j being in range. I do carry out this check following years of experience with both systems that do and systems that don't do this check. I would hope that the routines I supply with this package will reduce your need to access elements of matrices so speed of access is not a high priority. Use iterators? -------------- Iterators are an alternative way of providing fast access to the elements of an array or matrix when they are to be accessed sequentially. They need to be customised for each type of matrix. I have not implemented iterators in this package, although some iterator like functions are used for some row and column functions. Memory management - reference counting or status variable? ---------------------------------------------------------- Consider the instruction X = A + B + C; To evaluate this a simple program will add A to B putting the total in a temporary T1. Then it will add T1 to C creating another temporary T2 which will be copied into X. T1 and T2 will sit around till the end of the block. It would be faster if the program recognised that T1 was temporary and stored the sum of T1 and C back into T1 instead of creating T2 and then avoided the final copy by just assigning the contents of T1 to X rather than copying. In this case there will be no temporaries requiring deletion. (More precisely there will be a header to be deleted but no contents). For an instruction like X = (A * B) + (C * D); we can't avoid one temporary being left over, so we would like this temporary deleted as quickly as possible. I provide the functionality for doing this by attaching a status variable to each matrix. This indicates if the matrix is temporary so that its memory is available for recycling or deleting. Any matrix operation checks the status variables of the matrices it is working with and recycles or deletes any temporary memory. An alternative or additional approach would be to use delayed copying. If a program requests a matrix to be copied, the copy is delayed until an instruction is executed which modifies the memory of either the original matrix or the copy. This saves the unnecessary copying in the previous examples. However, it does not provide the additional functionality of my approach. It would be possible to have both delayed copy and tagging temporaries, but this seemed an unnecessary complexity. In particular delayed copy mechanisms seem to require two calls to the heap manager, using extra time and making building a memory compacting mechanism more difficult. Evaluation of expressions - use two stage method? ------------------------------------------------- Consider the instruction X = B - X; The simple program will subtract X from B, store the result in a temporary T1 and copy T1 into X. It would be faster if the program recognised that the result could be stored directly into X. This would happen automatically if the program could look at the instruction first and mark X as temporary. C programmers would expect to avoid the same problem with X = X - B; by using an operator -= (which I haven't provided, yet) X -= B; However this is an unnatural notation for non C users and it is much nicer to write X = X - B; and know that the program will carry out the simplification. Another example where this intelligent analysis of an instruction is helpful is in X = A.i() * B; where i() denotes inverse. Numerical analysts know it is inefficient to evaluate this expression by carrying out the inverse operation and then the multiply. Yet it is a convenient way of writing the instruction. It would be helpful if the program recognised this expression and carried out the more appropriate approach. To carry out this "intelligent" analysis of an instruction matrix expressions are evaluated in two stages. In the the first stage a tree representation of the expression is formed. For example (A+B)*C is represented by a tree * / \ + C / \ A B Rather than adding A and B the + operator yields an object of a class "AddedMatrix" which is just a pair of pointers to A and B. Then the * operator yields a "MultipliedMatrix" which is a pair of pointers to the "AddedMatrix" and C. The tree is examined for any simplifications and then evaluated recursively. Further possibilities not yet included are to recognise A.t()*A and A.t()+A as symmetric or to improve the efficiency of evaluation of expressions like A+B+C, A*B*C, A*B.t() [t() denotes transpose]. One of the disadvantages of the two-stage approach is that the types of matrix expressions are determined at run-time. So the compiler will not detect errors of the type Matrix M; DiagonalMatrix D; ....; D = M; We don't allow conversions using = when information would be lost. Such errors will be detected when the statement is executed. How to overcome an explosion in number of operations ---------------------------------------------------- The package attempts to solve the problem of the large number of versions of the binary operations required when one has a variety of types. With n types of matrices the binary operations will each require n-squared separate algorithms. Some reduction in the number may be possible by carrying out conversions. However the situation rapidly becomes impossible with more than 4 or 5 types. Doug Lea told me that it was possible to avoid this problem. I don't know what his solution is. Here's mine. Each matrix type includes routines for extracting individual rows or columns. I assume a row or column consists of a sequence of zeros, a sequence of stored values and then another sequence of zeros. Only a single algorithm is then required for each binary operation. The rows can be located very quickly since most of the matrices are stored row by row. Columns must be copied and so the access is somewhat slower. As far as possible my algorithms access the matrices by row. An alternative approach of using iterators will be slower since the iterators will involve a virtual function access at each step. In fact, I provide several algorithms for operations like + . If one is adding two matrices of the same type then there is no need to access the individual rows or columns and a faster general algorithm is appropriate. Generally the method works well. However symmetric matrices are not always handled very efficiently (yet) since complete rows are not stored explicitly. The original version of the package did not use this access by row or column method and provided the multitude of algorithms for the combination of different matrix types. The code file length turned out to be just a little longer than the present one when providing the same facilities with 5 distinct types of matrices. It would have been very difficult to increase the number of matrix types in the original version. Apparently 4 to 5 types is about the break even point for switching to the approach adopted in the present package. Using const ----------- The memory management scheme introduces a problem when a matrix is declared const. Because an operator may want to recycle the memory of its operands these operands cannot be declared const. It isn't reasonable for a temporary matrix to be declared const. However, I don't know how to tell this to the C++ compiler. One possibility is to provide alternative versions of the operators for operands declared const. But then one gets the explosion in the number of operators. My solution is to include versions of the initialisers for matrices declared const. Otherwise, you need to use A.c() in place of A if A is declared const and you wish to use it in an expression. A calculus of matrix types -------------------------- The program needs to be able to work out the class of the result of a matrix expression. This is to check that a conversion is legal or to determine the class of a temporary. To assist with this, a class MatrixType is defined. Operators +, -, *, >= are defined to calculate the types of the results of expressions or to check that conversions are legal. Error handling -------------- The package does not have graceful exit from errors. All errors are treated as fatal. Originally I thought I would wait until exceptions became available in C++. This now seems to have been delayed. In any case I don't think exceptions will solve all the problems. Some clean up of objects on the heap will often be required before one can exit via an exception. There are four categories of errors: Programming error - eg illegal conversion of a matrix, subscript out of bounds, matrix dimensions don't match; Illegal data error - eg Cholesky of a non-positive definite matrix; Out of space error - "new" returns a null pointer; Convergence failure - an iterative operation fails to converge. For the first two of these, it might be sensible to terminate a program. For the second two, one does want to return control to the user's program in a convenient manner. I don't know a good way of doing this, especially before exceptions are implemented. Band and sparse matrices ------------------------ The package does not yet support band or sparse types. At present the package assumes that the structure of a matrix is determined by its class and dimensions. This is not sufficient for band and sparse matrices. For band matrices one also needs to know the upper and lower band widths. For sparse matrices there is going to be some kind of structure vector. These are going to have to be calculated for the results of expressions in much the same way that types are calculated. In addition, a whole new set of row and column operations would have to be written for sparse matrices. However the present ones will be fine for band matrices. Band and sparse matrices are important for people solving large sets of differential equations. Sparse matrices are also important for statistical and operational research applications. ------------------------------------------------------------------------------- Matrix package problem report form ---------------------------------- Version: ...............newmat03 Date of release: .......Nov 28th, 1991 Primary site: ..........comp.sources.misc Downloaded from: ....... Your email address: .... Today's date: .......... Your machine: .......... Compiler & version: .... Describe the problem - attach examples if possible: Email to robert@am.dsir.govt.nz or Compuserve 72777,656 -------------------------------------------------------------------------------